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Abstract 



p 

We describe a novel approach to statistical learning from particles tracked while 
NO ', moving in a random environment. The problem consists in inferring properties 

of the environment from recorded snapshots. We consider here the case of a fluid 
seeded with identical passive particles that diffuse and are advected by a flow. Our 
approach rests on efficient algorithms to estimate the weighted number of possi- 
ble matchings among particles in two consecutive snapshots, the partition func- 
q ' tion of the underlying graphical model. The partition function is then maximized 

over the model parameters, namely diffusivity and velocity gradient. A Belief 
Propagation (BP) scheme is the backbone of our algorithm, providing accurate 
results for the flow parameters we want to learn. The BP estimate is addition- 
ally improved by incorporating Loop Series (LS) contributions. For the weighted 
matching problem, LS is compactly expressed as a Cauchy integral, accurately 
estimated by a saddle point approximation. Numerical experiments show that the 
quality of our improved BP algorithm is comparable to the one of a fully poly- 
nomial randomized approximation scheme, based on the Markov Chain Monte 
Carlo (MCMC) method, while the BP-based scheme is substantially faster than 
the MCMC scheme. 



1 Introduction 



$_h ' Graphical model approaches to statistical learning and inference are widespread in many fields of 

science, ranging from machine learning to bioinformatics, statistical physics and error-correction. 
Such applications often require evaluation of a weighted sum over an exponentially large number of 
configurations — a formidable #P-hard problem in the majority of cases. 

In this paper we focus on one such difficult problem, which occurs when tracking identical particles 
moving in a random environment. As long as particles are sufficiently dilute, their tracking in two 
consecutive frames is rather straightforward. When the density of particles and/or the acquisition 
time increase, many possible sets of trajectories become statistically compatible with the acquired 
data and multiple matchings of the particles in two consecutive snapshots are likely. Despite of these 
uncertainties, one expects that reliable estimates of the properties of the environment should still be 
possible if the number N of tracked particles is sufficiently large. This is the problem that we want 
to address here. 

The nature of the moving particles and their environment are not subject to particular restrictions, 
e.g. they might move actively, such as living organisms, or passively. Here, we shall consider the 
case of a fluid seeded with passive particles, a problem arising in the context of fluid mechanics 
experiments. Given a statistical model of the fluid flow with unknown parameters, along with the 
positions of N indistinguishable particles in two subsequent snapshots, one aims at predicting the 
most probable values of the model parameters. This task is formally stated in Section 2 as searching 
for the maximum of a weighted sum over all possible matchings between particles in the two snap- 
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shots. The problem turns out to be equivalent to computing the permanent of a non-negative matrix, 
known to be a #P-complete problem [11]. The main contribution of this paper is an efficient and 
accurate algorithm of Belief Propagation (BP) type for calculating the permanent for the class of 
weight matrices arising from the particle tracking problem. The BP algorithm seeks a minimum of 
the Bethe Free Energy [13] for a suitable graphical model. The graphical model is a fully connected 
bipartite graph: nodes are associated with the measured particles, edges are weighted according to 
the model of the flow transporting the particles and constraints enforce the condition that exactly 
one edge per node is active. It is known that BP gives the exact result for the maximum likelihood 
version of the problem (finding a maximum weight matching) in spite of multiple loops character- 
izing the graphical model [2]. The BP algorithm for the matching problem is derived and discussed 
in Section 3. 

BP equations could be understood as a re-parametrization, or gauge transformation, of factor func- 
tions in the graphical model [12]. Furthermore, BP solutions also provide an explicit representation 
of the exact partition function in terms of the so-called Loop Series [5, 6]. Our main technical result 
is the derivation of a compact expression and efficient approximation for the Loop Series in the prob- 
lem of weighted particle matching. This is done in Section 4, where the Loop Series is expressed 
in terms of an 2A-th order mixed derivative of an explicit functional, reduced to 2A-dimensional 
Cauchy integral and finally estimated by a saddle-point approximation. Section 5 describes empir- 
ical results demonstrating the performance of bare BP and the saddle-point improved BP in com- 
parison with a (simplified) fully polynomial randomized approximation scheme for computing the 
permanent [8]. Our improved BP achieves comparable accuracy, with significant gains in terms of 
speed. As the number of particles tracked in experiments is typically large (order tens of thousands) 
we argue that our approach is both useful and promising for applications. 



2 Particle tracking problem 



An important part of modern experiments in fluid mechanics is based on tracking of pre-seeded 
particles by sophisticated optical methods [1]. If particles are sufficiently small and chosen of ap- 
propriate (mass) density, their effect on the flow is essentially negligible and one can safely assume 
that they are passively transported by the flow. The (number) density of particles is usually rather 
high and a single snapshot typically contains a large number of them. The reason is that the smallest 
scales of the flow, which is generally turbulent, ought to be resolved. Two decades in a three dimen- 
sional flow require to follow at least one million, 10 2x3 , particles. Furthermore, turbulence is quite 
effective in rapidly transporting particles so that the acquisition time between consecutive snapshots 
should be kept small. Modern cameras have impressive resolutions, in the order of tens of thousands 
frames per second, yet the flow of information is huge: ~ Gigabit /s to monitor a two-dimensional 
slice of a (10cm) 3 experimental cell with a pixel size of 0.1mm and exposition time of 1ms. This 
extremely high rate makes it impossible to process data on the fly, unless very efficient algorithms 
are developed. 

Previous points motivate the development of a novel set of algorithmic tools for fast and efficient 
particle tracking. One key element is incorporating statistical models of the environment where 
particles are transported and tracked. For turbulent flows, modeling proceeds as follows. Consider 
N particles from the same time frame, labeled by i = 1, • • • , N and positioned at the set of points 
Xi, such that the typical distance between neighboring particles is smaller then the viscous scale 
of the flow. Then, Lagrangian particles evolve according to the set of stochastic equations, pi = 
U + Spi + £i, where pi are particle displacements on a line (generalization to multiple dimensions 
is straightforward) measured with respect to a reference point; U and S are the large-scale mean and 
gradient of the velocity field; £j (t) is the stochastic zero-mean Gaussian Langevin noise, describing 
molecular diffusivity, defined by its correlation function: (£i(ii)£j(i2)} = nSijS(ti — £2)- Particles 
are indistinguishable and the matching problem consists in assigning each particle from the original 
frame Xi = pi(Q) to particles in the subsequent frame y % — p»(A). Even if the flow parameters, 
U and S, and the diffusion coefficient, k, were known and frozen in time (the latter is a reasonable 
assumption provided the acquisition time A is sufficiently small), the matching cannot be identified 
with absolute certainty due to the stochastic nature of diffusion. The problem can be statistically 
modeled considering all possible particle matchings a between two successive frames and weighting 
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them according to 

\ «(cxp(2S)-l) J 

Here, of G {0, 1} is a Boolean variable indicating absence/presence of matching between x, and 
yJ, the vector a = (oi\i,j = ,7V) and = J]j :<KEj °i 1) Hi s (T,j <4> x ) enforces the 
constraints for a perfect matching (all particles match with exactly one particle in the other frame). 
For simplicity, U = (the drift common to all particles is subtracted) and time is rescaled to have 
A = 1. The partition function Z is the weighted sum over all possible matchings and p(a)/Z is 
their normalized probability distribution. By construction, the partition function is the permanent of 
the N x N positive matrix, p = (pj \i, j = 1, • • • , N), i.e. Z = per(p). 

Our goals are: (1) For given parameters S, K and the set of particle positions x and y in two subse- 
quent frames, have an algorithm for finding (a) the most probable matching, (b) marginal matching 
probabilities for any two particles from different frames, which is equivalent to computing the par- 
tition function. (2) Learn and provide reliable estimates of the model parameters S, k. 

Problem (la) is solved by the auction exact polynomial algorithm [3]. Conversely, problems (lb) 
and (2) belong to the #P-complete class, i.e. are likely to be exponentially complex, and we then 
aim at developing an efficient and systematically improvable heuristics. Our approach is based on 
the observation, made in [2], that a BP scheme equivalent to the auction algorithm can be formulated 
for (la), in spite of the underlying fully connected bi-partite graph with multiple loops (see also [4]). 
Notice that problem (la) is the Maximum-Likelihood version of (lb). We solve the problem (2) by 
taking the best possible estimate for Z at given values of S and k, and then maximizing the result 
over these parameters. We observed empirically that estimates based on Expectation Maximization 
(EM) algorithm [7] do not ensure accurate learning of the flow parameters S and k in some of 
the scenarios of interest, in particular the one with diffusion only. Conversely, BP gives accurate 
results in terms of the position of the maximum w.r.t the parameters, although the estimate for 
Z is often orders of magnitudes wrong. To further improve on this, we apply, in Section 4, the 
general BP-based Loop Calculus approach developed in [5, 6] to the perfect matching problem. 
This significantly improves the estimates of Z, especially in difficult cases when uncertainties in the 
matchings are significant. 

To the best of our knowledge, particle tracking as a learning problem - not to mention the algo- 
rithmic developments based on contemporary inference methods presented below - is novel and it 
was not discussed previously (see [9] for a survey of algorithms currently used in fluid mechanics 
experiments). 



3 Belief propagation and Bethe free energy 

For a model with states a having weight p(a) (as in (1)), the convex functional 

W?)} = In ^ (2) 

a 

has a single minimum, at b(a) = p{a)/Z (under the normalization condition ^2 s b(a) = 1), and 
the corresponding value of the functional T is the free energy, — In Z. As shown in [13], the Bethe 
free energy approximation and BP equations stem from (2) by considering an ansatz of the form 

Ma) w ^ — ; . (3) 

n w) 6?(W) 



Vectors <?; = {aj \j = 1, • • • , N} and cr- 7 = {aj \i = 1, • • • , N} are allowed to take any of the N 
possible values (0, • ■ • , 0, 1, 0, • • ■ ,0) with exactly one nonzero entry. Beliefs b\ (of), bi{5i), V{a^) 
satisfy for any i and j the consistency relation for marginal probabilities: bl(af) = bi(<Ji) = 

Ylg3\o-i where J2g .u denotes the sum over all possible values of the vector di keeping 
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fixed the value of the component a\. Eq. (3) is exact for a tree and serves as an approximation for 
graphs with loops, e.g. for the fully connected bi-partite graph of our matching problem. Beliefs, as 
approximations for probabilities, should also satisfy the normalization conditions: V(i, j), 6^(1) + 
6^(0) = 1. Using the normalization and consistency conditions, we can express all beliefs via 
ft = tf, (1) an d obtain for the Bethe free energy and the normalization conditions 

Fb P {0) = E [ft ln 4 - C 1 - ft) Ml -ft)), (4) 



Pi 



Vz: £#'=1; Vj: >J K L. (5) 



A simple argument shows that constrained minima of (4) are either a perfect matching (beliefs are 
all zeros and N of them are unity) or they are attained in the interior of the domain. The latter is 
the case generally encountered in the situations of interest to us, i.e. where no statistically dominant 
matching is present. Minima in the interior are stationary points of J- bp, corresponding to the 
following set of equations: 

V(*,j): ^;(! n -,,'). (6) 

where the 2N Lagrangian multipliers fj, (chemical potentials) are determined by Eqs. (5). Note that 
the Bethe free energy is not convex, and thus multiple minima in the interior of the domain might be 
possible. Empirically, we never found more than one though. In the limit where only the Maximum 
Likelihood configuration is of interest, entropy terms are discarded and (4) reduces to Linear Pro- 
gramming, yielding optimal integer solution in accordance with [2, 4]. Another relevant remark is 
that convexity is restored for a modified expression of the free energy, where the minus sign of the 
second term in Eq. (4) is reversed. The latter expression follows from an integral representation for 
Z, approximated in a saddle-point way. This approximation overestimates the diffusion coefficient 
and we use its unique solution (easy to find numerically) as initial condition to the following iterative 
version of Eqs. (5,6): 

V(i,j) : ft(n+l) = Xft(n) + -. = ( l ~ X )p1 ; (7) 

p' i + {Y i kftM/l + Y.kPKn)/2-ft{n)y/{u i {n)vj{n)) 



V*: ,(» + l) = ^f, Vi: ^1) ^ , 



where the arguments of the /3's indicate the order of the iterations, Ui — cxp(/ii) and = cxp(/i J ). 
The damping parameter A (typically chosen 0.4-^0.5) helps with convergence. To ensure appropriate 
accuracy for solutions with /3's close to zero or unity we also insert a normalization step after Eqs. (7) 
but prior to Eqs. (8), making the following two transformations consequently, (a) \f{i,j): (3f — > 

ft I Efc ft' anc ^ (b) ^(*' j) : ft ~* ft I Efc ft- Numerical experiments show that this procedure 
converges to a stationary point of the Bethe free energy (4). 



4 Loop series, Cauchy integral and saddle-point approximation 

As shown in [5, 6], the exact partition function of a generic graphical model can be expressed in 
terms of a Loop Series (LS), where each term of the series is computed explicitly using the BP 
solution. Adapting this general result to the matching problem, bulky yet straightforward algebra 
leads to the following exact expression for the partition function Z defined in Eq. (1): 

Z = Z BV *z, z^l+^rc, r c = (n(l-%)) f IIC 1 -^)] II 7^7 • < 9 > 

Here, the Bethe free energy Tbv = — \nZg-p, the variables (3 are in accordance with Eqs. (5,6), 
and C stands for an arbitrary generalized loop, defined as a subgraph of the fully connected bi- 
partite graph with all its vertexes having degree of connectivity > 1. The q t (or q J ) in Eq. (9) is 
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the C-dependent degree of connectivity of nodes, i.e. 5, = ^ 



{3\(i,l)dC} 



1 and q J = ^ 



1. 



According to Eq. (9), loops with even/odd number of vertexes give positive/negative contributions 
rc- Therefore, the series is not positive definite, which is also consistent with the fact that Zg-p in 
general does not provide a lower bound for the exact partition function. (In some special cases, e.g. 
for the model studied in [10], all terms in the series are known to be positive and thus Zpp < Z.) In 
all cases of the weighted matching problem we have experimented with, we have empirically found 
that the inequality Zbp < Z still holds. Let us finally notice an important special feature of the 
weighted matching problem: for any generalized loop C, its individual contribution \rc\ < 1- The 
proof can be found in the Appendix. 

Eq. (9) allows for the following compact representation in terms of 2 N-th order mixed local deriva- 
tive of an explicit function of 2N variables 



(10) 




(id 



where p = (pi, • • • , Pni P > ' • ' j P ) are auxiliary variables. However, calculating the 2A r -order 
mixed derivative exactly is a task of exponential complexity and one wonders whether the mixed 
derivative can be approximated efficiently. Partial answer to this question is given below. 

Using the Cauchy integral representation for the first-order derivative of an analytic function (Z 
is analytic over pi, p^ with finite real parts), Eqs. (10,11) can be recast as the following contour 
integral: 



(f> exp (- 



■em 



II dp II '//' 



(2tti) 



2^ 



g(p)= 21npi + ^bV'-lnZ, (12) 



where T p is a direct product of 27V close contours circling clockwise the origin p = 0, where 
derivatives in (10) are to be computed. We observe that each integral over an individual p variable in 
Eq. (12), say pi, has a pole at pi = and essential singularities at pi = ±00. Notice also that 9(p) 
is a concave function of p in each one of the 2 2N quadrants, i.e. where components of p are finite, 
real and have a definite sign. Therefore, it is natural to shift the contour of integration to one of 2 2N 
maxima of Q{p), 



Vi 




Vj: 7 = 1 ~^ 



Pi 



(13) 



and orient the contour along the direction of steepest descent from the saddle point. Once the sign 
of each component of p is fixed, approximating the respective solution of Eqs. (13) numerically is 
straightforward due to the concavity of Q. We shall enumerate the various maxima by the index s. 

The saddle-point approximation of Eq. (12), accounting for Gaussian integral corrections about all 
the maxima of Eq. (13), yields 



cxp(-e(^ ) )) _ 



2^™ 

E 

s =i (27r) JV A/det(A( s )) 



]Texp(-^>(e< 



8) 



(14) 



where A' s ' is the Hessian of Q{p) at the saddle-points. Corrections to each term in Eq. (14) 
are measured in terms of higher-order terms, the leading (fourth-order) being estimated as 
Qi S) = -iL^^^AWj^^ttAW)- 1 )^, where is the tensor of fourth- 

order derivatives of Q(p). The improved approximation for the partition function becomes z « 
s exp(— Gspitf^) — Giio^^))- One reason to account for the fourth-order term is that the ra- 
tio 1 1?4 / C? sp I gives a standard measure of the saddle -point validity. In cases where the saddle- 
approximation becomes asymptotically exact, the ratio approaches zero. A weaker condition, 
I ^4 / S p I < L suffices for a heuristically reasonable approximation. 
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Figure 1: Numerical results comparing BP, Loops series improvement and MCMC. Left panel: 
diffusion only; right panel: diffusion and advection. Red curve are BP results, blue and green are 
Loop series improvements, and black curve is MCMC. 

We also expect that the typical order of magnitude of the 2 2N terms Qs'p and Q± grows as O(N). 
Therefore, the sum in s over the saddle-points might be dominated by a single term in the limit of 
large N. As discussed in the next Section, we find empirically that such dominant term indeed exists 
and happens to correspond to the vector having all its components positive. Therefore, the sum 
over all the maxima, indexed by s, can be simplified by keeping only the dominant contribution 
correspondent to g 1 -^ . The expression that we have employed in numerical experiments discussed 
in the next Section reads finally, z w exp (— G sp ((f + ') — Qiiff^)) ■ 

5 Numerical results 

In this Section we compare the accuracy of BP based approximations for the partition function, Z, of 
the model (1) with MCMC simulations. As briefly stated in Section 2, computing Z for the weighted 
matching problem is equivalent to computing a permanent of a non-negative matrix, for which a 
Fully Polynomial Randomized Approximation Scheme (FPRAS) exists based on the Markov Chain 
Monte Carlo method [8]. We implemented the basic idea of FPRAS, with some simplifications 
applicable to our problem. This algorithm was used to assess accuracy of our approximations, but is 
orders of magnitude slower than the BP based approaches, and thus not applicable to large particle 
tracking problems. 

To study dependence of Z on n and S (see discussion of Section 2), we estimate Z at different values 
of K and S and compare the curves, searching for the maximum with respect to these parameters. 
Our BP simulations consist of the following steps for each value of k or S. First, we find solution 
of BP equations running the numerical scheme described in Eqs. (7,8) and calculating the resulting 
J- bp according to Eq. (4) . Second, we find the solution of the saddle-point Eqs. (13) in the various 
quadrants, calculate the covariance matrix A( s ) for this saddle solution and thus estimate the leading 
saddle-point correction Q sp in accordance with Eq. (14). Finally, we estimate the respective fourth- 
order correction, Q± . 

We compare respective contributions to the partition function associated with saddle-points with 
different choices of the signs. In all experiments where the typical overlap among particles is sig- 
nificantly smaller than the total number of particles we find that the contribution with all signs + 
dominates. Moreover, the gap separating the leading + contribution and other contributions is sig- 
nificant and grows linearly with N. This allows us to ignore all other saddle-point contributions 
but the + one. (In cases of moderate N we have made the exhaustive comparison. In general, we 
compared the all + contribution with all — contribution, contributions with a limited number of 
signs flipped and also with signs generated randomly.) We observe that the saddle-point validity 
conditions holds reasonably well if N is sufficiently large, at N = 100 the ratio \Gi/Q S p\ is typi- 
cally 0.1 -r 0.4. Since \Qi/Q sp \ — ► when S — > +oo, we find that the saddle-point is very accurate 
(possibly asymptotically exact) in this limit. 
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Results of numerical simulations for number of particles TV = 100 are shown in Figure 1. The plots 
show results for one set of particle positions each, but we found that variability between scenarios 
decreases with increasing number of particles, and one scenario is thus sufficient to characterize the 
main trends. 

The left panel shows results for a situation with diffusion only, with actual diffusivity n ac t = 1.0 
used for generating the particle positions. The x-axis spans values of the governing parameter k 
(pretended not to be known), and the y-axis shows configuration weight per particle (what is shown 
is In Z/N). The black curve are results of the MCMC simulation, which are close to the true values 
of Z. The curve indeed peaks around the correct value n = K ac t = 1.0, although the maximum is 
rather shallow (the lower K ac j the more pronounced is the maximum). The red curve corresponds 
to results obtained by BP only, and shows that BP severely underestimates the partition function, 
although its maximum is at the right value. The remaining blue and green curves are the two saddle 
corrections discussed in Section 4 As can be seen, the saddle corrections significantly improve 
the BP estimate. The time to compute each point in the plot using the BP based scheme is w 5sec, 
while the MCMC algorithm takes pa lOmins. 

The right panel of Figure 1 shows results for a scenario with both advection and diffusion. The 
x-axis now spans values of the velocity gradient S, and y-axis is again In Z/N. The actual velocity 
gradient used for generating the particles was 5 ac t = —1.0. The four curves show similar main 
trends: BP underestimates Z and loop corrections provide very tight belt around the MCMC results. 
The main difference is that in the case when advection is present, the peak is very well pronounced 
and all methods give extremely accurate answer for the velocity gradient S. In this case, the running 
times were again w 5sec for the BP scheme, but w 5hours for MCMC per point. 

6 Conclusions 

We have presented new computational tools for particle tracking, based on Belief Propagation and 
Loop Series and compared them to the Markov Chain Monte Carlo scheme for the estimation of 
the permanent [8]. We have specifically considered tracking of passive particles in fluid dynamics 
experiments. The methods are quite general though and applications to other tracking problems, 
e.g. to self-propelling biological objects, are possible and will be pursued. Our long-term goal is 
to develop computational tools effective enough to permit on-the-fly processing of particle tracking 
frames. Algorithms presented here ensure an excellent accuracy and the BP-improved scheme is 
already orders of magnitudes faster than the MCMC scheme. 

Appendix 

Here we prove the following important property of the loop series for perfect matching: 

Proposition 1. In the loop calculus expansion of the graphical model for perfect matching prob- 
lem (9), \rc\ < 1 for all generalized loops C. 

Proof. We rewrite expression (9) as: rc = (Iliec V'ijC?) flX/eC ^CJ> wnere V^c = (1 — 

1i) Y\.{j\(i j)ec} V ^M/0-~ Pi)' an ^ analogously for ip J c . We will use the fact that in a fixed point 
of our BP scheme J^i^i = Si Pi = !■ We proceed by showing that \ipi-,c\ is maximized for 
f3 3 > = ... = = l/q u and that even for such beliefs \ipi-c\ < 1- The situation is analogous for 
x/jq. It then follows that \rc\ < 1 as desired. 

'We find that the quality of saddle-approximation decreases with k, when the number of "polarized" match- 
ings, i.e. those with j3j — » 1, increases. Polarized beliefs, correspondent to almost committed matchings, do not 
contribute significantly to the Loop Series yet their respective contributions are misrepresented in the saddle- 
approximation. To compensate for this caveat of the saddle-approximation we decimated the original graph, 
reducing it to a subgraph with all particles involved in "polarized" matchings (and all edges associated with 
them) pruned out. For the numerical experiments in Figure 1, we used the polarization criterion, j3\ > 0.01. 
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Let us fix a generalized loop C and a node i with degree % > 2 in C and seek values of (3? 1 = 
. . . = S [0, 1] s.t. Y?k=i Pi° ^ 1 that maximize \ipi-c\- First, by observing that (3/(1 - (3) 
is a growing function of /?, we see that indeed Y^k=i Pi k = 1 must be the case. Focusing only on 
the subexpression of with (3 in it, we have the following constrained optimization problem: 

max/3 Ylj-(i j)ec { s - t - Sfe=i /^i* = ^ We a PPty me method of Lagrangian Multipliers, and 
arrive at the following conditions for extrema: 

Taking any two j\ ^ j2, we derive p 32 (l — (3f 2 ) = (3^(1 — from the above. This, being a 
quadratic equation in ffl x , has exactly two solutions, and it is not difficult so see that they must be 
(3? 1 = (3 32 or Pf 1 = 1 — (3 32 . These are the conditions on /3s for an extremum of the product. 

To incorporate the (1 — qt) subexpression of tjj^c, let us deal with the special case when q L = 2 
separately: in this case, since /3f 1 + (3f 2 — 1, we have that the whole product we maximize is equal 
to 1, and thus \ipi\c\ = 1. For g, > 2: if, for any pair Pf 1 = 1 — (3 d2 , then all other (3s 

(other than Pf 1 and (3j 2 ) must be due to J2t=i P\ k = ^- This means that the product is zero, and 
therefore ipi-.c = 0. If, on the other hand, fjP- = (3 J2 for all pairs ji, j 2 , we have that (3j = 1/qt- 
In this situation, the question is whether — 1) 1_ 9 ; / 2 < 1, which is true for any > 2. In fact, 
we have shown that the expression (gj — l) 1-1 ?*/ 2 is an upper bound on \ipi-c\ f° r an y node i in any 
generalized loop C. 

□ 
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